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Abstract 

We will present a complete set of equations, in the form of an Einstein- 
Bianchi system, that describe the evolution of generic smooth lattices 
in spacetime. All 20 independent Riemann curvatures will be evolved 
in parallel with the leg-lengths of the lattice. We will show that the 
evolution equations for the curvatures forms a hyperbolic system and 
that the associated constraints are preserved. This work is a gener- 
alisation of our previous paper [1] on the Einstein-Bianchi system for 
the Schwarzschild spacetime to general 3+1 vacuum spacetimes. 



1 Introduction 



In a series of papers we have shown that the smooth lattice method works 
remarkably well for simple spacetimes such as the Schwarzschild spacetime in 
various slicings [1, 2], the maximally sliced Oppenheimer-Snyder spacetime 
[3], the vacuum Kasner cosmologies [1] and for constructing Schwarzschild 
initial data [5]. The equations are simple and require little computational 
sophistication to achieve stable and accurate results. The real test of the 
method however must be in the context of generic spacetimes. This paper is 
a first step in that direction. 
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The logic behind the smooth lattice approach is quite simple. We assume 
that we are given a lattice, built from a large collection of interconnected 
vertices, and where each path that connects a pair of vertices is taken to be 
a geodesic segment of the spacetime. The only data that we are given for 
the lattice is the connection matrix (which describes the topology as a list of 
pairs of connected vertices) and the lengths of each geodesic segment (which 
describes the metric properties). In this picture we are assuming that the 
lattice geometry is a close approximation to some underlying smooth geome- 
try. The question that (should) spring to mind is - Given the leg lengths on a 
lattice, how do we compute the Riemann curvatures? We will return to this 
important question in just a moment, but for now let us suppose we have a 
suitable algorithm by which we can accurately compute the Riemann curva- 
tures. It is then a simple matter to impose the vacuum Einstein equations* 
which in turn will impose constraints^ on the leg-lengths. This furnishes us 
with a discrete set of equations for the leg-lengths. Solving these equations 
will yield a discrete solution of the vacuum Einstein equations. 

We now return to the question of how to recover the Riemann curvatures 
given the set of leg-lengths. In one of our earlier papers [5] we argued that if 
the lattice was sufficiently well refined then a local Riemann normal coordi- 
nate frame could be constructed in the neighbourhood of any vertex extend- 
ing to include, at least, the immediate neighbouring vertices. We called this 
neighbourhood the computational cell for the vertex (for lattices built from 
tetrahedra this would consist of the tetrahedra attached to the vertex). In 
this computational cell we can expand the metric as a power series [( i] around 
the central vertex 



where L is a typical length scale for the computational cell. The requirement 
that the legs are geodesic segments leads, after some detailed calculations 
[ti], to the following equation 



where Aa;^- := — . The approach advocated in [ ] was to use this 
equation to extract the Riemann curvatures from the lattice. This may sound 
simple but there number of troubling issues. 

*For pure pedagogy we will restrict the discussion to vacuum spacetimes. 
'''Not to be confused with any constraints that may exist at the continuum level, for 
example the ADM constraints. 
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The first issue concerns the coordinates. How do we compute coordinates for 
each vertex? Some can be set by simple gauge transformations (e.g., the ori- 
gin can be tied to the central vertex) while the remainder must be computed 
from the lattice data (i.e., the leg-lengths). This forces us to view the above 
equations (1.2) as a coupled system for the curvatures and the coordinates. 

The second issue is one of accountancy - do we have enough equations to 
compute the curvatures and the coordinates? For most lattices (in 3 and 
higher dimensions) the legs out number the coordinates and curvatures 
R^lau|3■ As an example, the computational cell used in our earher paper [5] 
contained 78 legs and 19 vertices. Thus we had 78 equations for 6 curvatures 
and 57 coordinates (of which 6 can be freely chosen). There are at least two 
ways to handle this over supply of information. We can either form linear 
combinations of the above equations (1.2) to produce a reduced system in 
which the number of equations matches the number of unknowns. Or we 
can include a sufficient number of higher order terms in the Taylor series 
so as to produce a consistent set of equations. This later approach has the 
possible benefit of producing higher order approximations for the R^avis but 
at considerable extra expense. In both instances we still have a large coupled 
non-linear system of equations to solve at each vertex and at each time step. 
This is a considerable computational challenge. 

Another important issue is one of uniqueness - how many distinct solutions 
can we find for the and -R^a^^? The equations are non-linear and thus it 
is conceivable that more than one solution could be found. Do the solutions 
form a continuous family or are there only a finite set of solutions? How 
would we choose between these solutions? In our earlier paper [5] we resolved 
these problems by extending the lattice data to include the angles between 
each pair of legs attached to the central vertex. This allowed us to obtain 
an explicit and unique solution for all of the coordinates in a computational 
cell. It also had the added bonus of decoupling the coordinates from the 
curvatures - we could calculate all of the coordinates before computing the 
curvatures. The price we paid for this improvement was a significant increase 
in the number of data to be evolved. Where previously we had 78 legs per 
computational cell, now we had a further 33 angles. 

However, there is a final issue which is much more serious than those just 
mentioned. To obtain O (L) accurate estimates for the curvatures, the coor- 
dinates must be computed to at least O {L^) accuracy (i.e., the errors must be 
no worse than O {L^))- This follows by inspection of equation (1.2). Suppose 
the error in is O {L"') for some a > 0. This error will couple with the first 
term on the right hand side of (1.2) to introduce an error of O {L"-^^). But 



3 



the curvature terms are O (L^) and will dominate the error term only when 
a > 4. Admittedly this is a somewhat naive analysis as it takes no account of 
the smoothness of the underlying geometry which might ensure that various 
lower order terms cancel (see for example the role smoothness plays in es- 
tablishing the truncation errors in centred finite-difference approximations). 
But in the absence of an explicit algorithm we are unable to demonstrate 
that such cancellations do occur-'-. The upshot is that if we persist with any 
of the variations suggested above we must design a solution strategy that 
guarantees, without invoking smoothness, that the errors in the coordinates 
are no worse than 0{L^). Despite our best efforts, we have not found a 
reliable solution to this problem. 

These issues are not altogether new nor surprising and have proved to be a 
niggling concern throughout the development of the smooth lattice method. 
The only working solution that we have found (there may be others) is to 
surrender some (or all) of the main equations (1.2) in favour of the Bianchi 
identities. In all of our papers [1, 2, 3, 4, 5] we used a combination of the 
Bianchi identities and the geodesic deviation equation in 1 -|- 1 spacetimes. 
The results were very encouraging. This was a hybrid scheme^ and we at- 
tributed its success to the introduction of the Bianchi identities. This is 
the motivation for the present paper - Can we use the Bianchi identities to 
compute all of the Riemann curvatures in a 3 + 1 spacetime? We should 
emphasise that there is one important difference between what we propose 
here and our previous work. In this paper we will use the full set of Bianchi 
identities to evolve all 20 independent Riemann curvatures. In contrast, in 
our 1 + 1 experiments we used one Bianchi identity to compute one spa- 
tial curvature (i.e., a purely 3-dimensional computation within one Cauchy 
surface). 

Why should we believe that this use of the Bianchi identities will overcome 
the issues described above? Simply, it allows us to use lower order approxi- 
mations for the vertex coordinates (even flat space approximations) without 
compromising the quality of the estimates of the curvatures. We will return 
to this point after we have presented the full set of evolution equations. 

-'-Though the introduction of angles does produce an exphcit algorithm its analysis is 
too un-wieldily to be of any use. 

^The geodesic deviation equation arises as a continuum limit of the smooth lattice 
equations [ j]. 
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2 Notation 



A typical computational cell will be denoted by Q. This will be a compact 
subset of the spacetime manifold. The central vertex of the cell will be 
denoted by O and the subset of Q obtained by the intersection of fl with 
the particular Cauchy surface that contains O will be denoted by u. We will 
describe u as the floor of Q. As Q has a finite extent there will be an image 
of u that defines the future end of Q. We will refer to this as the roof of Q. 
We will have little to reason to refer to the past end of Q but calling it the 
basement seems consistent. 

We will assume throughout this paper that the vertex world lines are normal 
to the Cauchy surfaces (i.e., zero drift, in the language of [1]). This may 
seem restrictive but in our experiments to date it has worked very well. 

Within Q we will employ two sets of vectors essential to the evolution of 
the lattice. The first set will be an orthonormal tetrad, denoted by e^, 
a = 1, 2, 3, 4, tied to the world line of O and aligned so that ei is the tangent 
vector to the world line of O. As we have assumed that the drift vector is 
everywhere zero this also ensures that ei is the future pointing unit normal to 
cu at O. Following convention, we will write n^^ as the unit normal to u though 
as just noted, this is identical to Ci. The second set of vectors will be based 
on the set of radial legs attached to O. Each leg will be of the form (oi) and 
we will use Vi to denote the vector that joins (o) to (i). Note that the Vi are 
neither unit nor orthogonal. Latin characters will always be used to denote 
tetrad indices while the spacetime indices will be denoted by Greek letters. 
Latin characters will also be used as vertex labels and where confusion might 
arise we will use subsets of the Latin alphabet with a,b,c, ■ ■ ■ h reserved 
for frame components while i,j,k,l,m will be reserved for vertex labels. 
Obviously this distinction will only be imposed for equations that contain 
both types of index. 

Each cell will carry a Riemann normal coordinate frame (an RNC frame), 
with coordinates x'^ = {t, x, y, z), tied to the central vertex and aligned with 
the tetrad. Note that this gives precedence to the tetrad over the coordinates. 
Coordinate components will be written as i?^,^ or for specific components as, 
for example, Rtx while for frame components we will use scripts characters 
TZab- The coordinates for a typical vertex (z) will often be written as x^ but 
on occasion we will have need to talk about the particular values for the 
in which case we will write {t,x,y,z)i or even xj, xf etc. 

Each RNC frame will be chosen so that at O the metric is diagonal, {g^u)o = 
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diag(— 1, 1, 1, 1). Both spacetime and tetrad indices will be raised and low- 
ered, at O, using the metric diag(— 1, 1, 1, 1). With these choices we see that 
the future pointing unit normal to the Cauchy surface at the central ver- 
tex O is just (n^)o = (1,0,0,0)^ while (n^)o = (-1,0,0,0). We also see 
that the tetrad Ca has components e^a = in this RNC frame. Note that 
e^ae^I' = 5a , e^ae^-u = S^^u, e^i = n'' and e^,^ = -n^. 



3 Evolving the leg- lengths 

The legs of the lattice are required to be short geodesic segments. Thus 
it should come as no surprise that the evolution of the leg-lengths can be 
obtained from the equations for the second variation of arc length. In an 
earlier paper [7] we showed that, for sufficiently short legs, these equations 
can be written as follows 



dt 



-2NK^,Axf^Ax';j + 0{L^) (3.1) 



d ( 1 d^, 



dt \N dt 



^ -= 2iV|,^A4Aa;f^. (3.2) 



For numerical purposes it is somewhat easier to rewrite these in the following 
form 

= (3.3) 

d P 

^ = -iV|,,Axr,Axg (3.4) 
- N {K^^K'^p - R.^^pn^'n'') Ax,^- Axg 

in which we have introduced the new variables Pjj, one per leg. The K^y 
can be obtained by a suitable weighted sum of equation (3.1) as described 
in section (7.1). We have also dropped the truncation terms as these are not 
used during a numerical integration. 

Clearly, the evolution of the leg lengths requires a knowledge of the Riemann 
curvatures and to that end we now present the evolution equations for those 
curvatures. 
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4 Evolving the Riemann curvatures. Pt. 1 



We know that there are only 20 algebraically independent Riemann curva- 
tures in 4 dimensions. So which should we choose? By a careful inspection 
of the algebraic symmetries of Rfiav^ we settled upon the following 

Rxyxy-i Rxyxz-i Rxyyzi Rxzxzi Rxzyzi Ryzyz 
R-txxyi R-txxzi R-tyxyi R-tyxzi R-tyyzi R-tzxyi R-tzyzi Rtzyz (^■■^) 
R-txtxt R-tytyt R-tztzt Rtxtyt Rtxtzt Rtytz 



4.1 Bianchi identities 

Our aim is to use the Bianchi identities to obtain evolution equations for the 
Riemann curvatures. We begin by writing down the Bianchi identities at the 
central vertex, where the connection vanishes, 

along with a contracted version of the same equation 

= g'^'^Rixav^a — Ral3,v + Rav,l3 (4.3) 

This pair of equations, along with the vacuum Einstein field equations, and 
a judicious choice of indices will provide us with all of the required evolution 
equations. This leads to the following 14 differential equations 






R'xyxy.t 


Rtyxy.x 


^" Rtxxy,y 


(4.4) 





RxyxZjt 


Rtzxy,x 


~l~ Rtxxy,z 


(4.5) 





Rxyyz,t 


Rtzxy,y 


~l~ Rtyxy,z 


(4.6) 





Rxzxz,t 


RtZXZjX 


~l~ Rtxxz,z 


(4.7) 





— Rxzyz,t 


Rtzxz,y 


~l~ Rtyxz,z 


(4.8) 





Ryzyz,t 


Rtzyz,y 


~l~ Rtyyz,z 


(4.9) 





Rtyxy,t 


Rxyxy,x 


~l~ Rxyyz,z 


(4.10) 





Rtxxy,t 


~l~ Rxyxy,y 


~l~ Rxyxz,z 


(4.11) 





— Rtzxy,t 


Rxyxz,x 


~ Rxyyz,y 


(4.12) 





Rtzxz,t 


RxZXZjX 


Rxzyz,y 


(4.13) 
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Rtxxz,t 


~l~ Rxyxz,y 


~l~ Rxzxz,z 


(4.14) 





Rtyxz,t 


Rxyxz,x 


~l~ Rxzyz, z 


(4.15) 





Rtzyz,t 


Rxzyz,x 


~ Ryzyz,y 


(4.16) 





Rtyyz,t 


Rxyyz,x 


~l~ Ryzyz,z 


(4.17) 



There are of course 20 independent Ruaup, 14 of which are subject to the 
above evolution equations while the remaining 6 can be obtained from the 
vacuum Einstein equations 






Rxx 


Rtxtx ~l~ Rxyxy ~l~ Rxzxz 


(4.18) 





= Ryy = 


Rtyty ~l~ Rxyxy ~l~ Ryzyz 


(4.19) 





Rzz 


Rtztz ~l~ Rxzxz ~l~ Ryzyz 


(4.20) 





Rxy 


Rtxty ~\~ Rxzyz 


(4.21) 





Rxz 


Rtxtz Rxyyz 


(4.22) 





Ryz 


Rtytz ~l~ Rxyxz 


(4.23) 



Though these are not differential equations they do, none the less, provide a 
means to evolve the 6 curvatures Rtxtx ^ Rtxty ■ ■ ■ Rtytz- 

The important point to note about this system of equations is that it is 
closed, there are 20 evolution equations for 20 curvatures. The source terms, 
such as Rxyxy, X, could be computed by importing data from the neighbouring 
cells, by an appropriate combination of rotations and boosts, and using a 
suitable finite difference approximation (see section (7) for more details)). In 
this way the lattice serves as a scaffold on which source terms such as these 
can be computed. 

4.2 Constraints 

In deriving the 20 evolution equations of the previous section we used only 6 
of the 10 vacuum Einstein equations. Thus the 4 remaining vacuum Einstein 
equations must be viewed as constraints. These equations are 






— Rtt — Rtxtx + Rtyty + Rtztz 


(4.24) 





Rtx Rtyxy ~\~ Rtzxz 


(4.25) 





Rty Rtxxy ~\~ Rtzyz 


(4.26) 





Rtz Rtxxz Rtyyz 


(4.27) 
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Finally, we have the following 6 constraints that arise from the Bianchi iden- 
tities. 




= R 



+ R., 



-R. 



(4.28) 
(4.29) 
(4.30) 
(4.31) 
(4.32) 
(4.33) 



So all up we have 20 evolution equations assembled from the 14 differential 
equations (4.4-4.17) and 6 algebraic equations (4.18-4.23) plus 10 constraints 
comprising 4 Einstein equations (4.24-4.27) and 6 Bianchi identities (4.28- 
4.33). This is a such a simple system that it allows simple questions to be 
explored and answered with ease. The questions that we will address are 

1. Are the constraints preserved by the evolution equations? 

2. Do the evolution equations constitute a hyperbolic system? 

For both questions the answer is yes and we shall now demonstrate that this 
is so. 

4.3 Constraint preservation 

In the following discussion we will assume that, by some means, we have 
constructed an initial data set for the 20 R^au/s- That is, the 20 R^au/s are 
chosen so that the 10 constraints (4.24-4.33) vanish at the central vertex of 
every computational cell in the lattice. 

We will also need the trivial result that 



which follows directly from equations (4.18,4.19,4.20,4.24). 

Consider now the constraint = Rtz- By assumption, this constraint is satis- 
fied on the initial slice. To demonstrate that it continues to hold throughout 
the evolution we need to show that = Rtz,t- From (4.27) this requires us to 
show that = Rtxxz,t + Rtyyz,t- Usiug (4.14,4.17) we see that 




(4.34) 




'xyyz,x 



-R. 
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however on the initial shce we also have, by assumption, (4.28) 

Rxyxy,z ~l~ Rxyyz,x Rxyxz,y 

which when combined with the previous equation leads to 

Rtxxz,t ~l~ Rtyyz,t {Rxyxy ~l~ Rxzxz ~l~ Ryzyz) ^ 

But by equation (4.34) we see that the right hand side is just —R^z/'^ and 
as i? = across the initial slice we also have that = R^z at every central 
vertex. This completes the proof. The two other constraints, = Rty and 
= Rtxi can be dealt with in a similar fashion. 

All that remains is to show that = Ru is conserved. We proceed in a 
manner similar to the above. First we use Ru = Rxyxy + Rxzxz + Ryzyz and 
then use equations (4.4,4.7,4.9) to compute the time derivative 

i^Rxyxy ~l~ Rxzxz ~l~ Ryzyz) i Rxyxy, t ~\~ Rxzxz, t ~\~ Ryzyz,t 

Rtyxy,x Rtxxy,y 
~l~ RtZXZjX Rtxxz,z 
~l~ Rtzyz,y Rtyyz,z 
Rtx,x ~l~ Rty,y ~l~ Rtz,z 

where the last line arose by inspection of equations (4.25-4.27). But = R^i, 
at every central vertex on the initial slice. Thus = Rfj,u,i, i = x,y,z on the 
central vertex which in turn shows that = Ru,t on the initial slice. 

A key element in the above proofs was the use of constraints based on the 
Bianchi identities. The question now must be - do the evolution equations 
preserve those constraints? The answer is yes which we will now demonstrate 
on a typical case. Consider the constraint (4.28) 

Rxyxy, z ~l~ Rxyyz,x Rxyxz,y 

We know this to be true on the initial slice and we need to show that the 
evolution equations (4.4-4.17) guarantee that it will be satisfied on all sub- 
sequent slices. The calculations follow a now familiar pattern, 

i^Rxyxy,z ~l~ Rxyyz,x Rxyxz,y) ,t Rxyxy ,tz ~\~ Rxyyz,tx Rxyxz,ty 

{,Rtyxy,x Rtxxy,y) ,z 
~l~ (^Rtzxy,y Rtyxy,z) ,x 
(^Rtzxy,x Rtxxy,z) ,y 

= 

The same analysis can be applied to the remaining constraint equations. 
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4.4 Hyperbolicity 



Our approach to proving hyperbolicity will be quite simple. We will manip- 
ulate the evolution equations (4.4-4.17) to demonstrate that each of our 20 
R^lau|3 satisfies the standard second order wave equation. 

Let us start with a simple example, equation (4.4). We take one further time 
derivative, commute the mixed partial derivatives and then use equations 
(4.10,4.11) to eliminate the single time derivative. This leads to 

Rxyxy ,tt Rxyxy,xx ^xyxy ;yy ~l~ Rxyyz,zx Rxyxz,zy 

However, we also have = Rxyxy,z + Rxyyz,x — Rxyxz,y, which allows us to 
reduce the last two terms of the previous equation to just —Rxyxy,zz- Thus 
we have 

Rxyxy ,tt Rxyxy,xx Rxyxy ,yy Rxyxy,zz 

This is the standard fiat space wave equation for Rxyxy A similar analysis 
shows that Rxzxz, Ryzyz, Rxyxz, Rxyyz dJid Rxzyz are also solutious of the wave 
equation. 

We now turn to the 8 R^avp in which the indices /iaz//3 contain just one t. 
The proof (that each such R^aup satisfies the wave equation) differs from the 
above only in the way the Bianchi identities are used. Applying the first few 
steps outlined above to equation (4.10) leads to 

Rtyxy,tt Rtyxy,xx Rtyxy,yy Rtyxy,zz 
~l~ Rtyxy,yy ~l~ Rtxxy,xy ~\~ Rtzxy,zy 

in which we have deliberately introduced the pair of terms Rtyxy,yy to aid in 
the following exposition. The last three terms can be dealt with as follows. 
First notice that 

Rtyxy,yy ~l~ Rtxxy,xy ~l~ Rtzxy,zy (^Rtyxy,y ~l~ Rtxxy,x ~l~ Rtzxy,z) ,y 

— {^R'^txy,fj.) ,y 

( Rty.x Rtx,y) ,y 

where in last line we have used the contracted Bianchi identity = R'^iyai3,fi — 
Rvi3,a + Rua,i3- But we kuow that = R^u at every central vertex, thus all 
of its partial derivatives will be zero and so the each term on the right hand 
vanishes leading to our desired result 

Rtyxy,tt Rtyxy,xx Rtyxy,yy Rtyxy,zz 
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Finally we note that the remaining 6 R^aup^ that is those that carry two t's in 
their indices, are linear combinations of the previous 14 R^aui3, see equations 
(4.18-4.23), and thus will also be solutions of the wave equation. Thus we 
have shown, as claimed, that all 20 R^aup satisfy the wave equation. 



5 Evolving the Riemann curvatures. Pt. 2 

There are two problems in the forgoing analysis. The first problem is that 
we chose a unit lapse function when presenting the evolution equations (4.4- 
4.17). We can easily remedy this problem by making a simple vertex depen- 
dent coordinate substitution t = Nt' in each of the evolution equations. 

The second problem is somewhat more of a challenge. It stems from the 
simple fact that each computational cell is local in both space and time and 
therefore no single RNC can be used to track the evolution for an extended 
period of time. We will have no choice but to jump periodically to a new RNC 
frame. But how might we do this? One approach goes as follows. Build, on 
the world line of a typical vertex, a pair of distinct but overlapping cells, with 
one cell lying slightly to the future of the other. Then evolve the curvatures 
in the frame of one cell into the overlap region followed by a coordinate 
transformation to import the newly evolved curvatures into the frame of the 
future cell. This completes one time step of the integration whereupon the 
whole process can be repeated any number of times along the vertex world 
line. A useful improvement on this is to use a local tetrad to construct 
scalars thus avoiding the need for explicit coordinate transformations when 
passing from one cell to the next. The price we pay for this is that we have 
to account for the evolution of the tetrad along the world line. As we shall 
see this is rather easy to do (essentially we project the tetrad onto the legs 
of the lattice). We will explore this method first on a simple example before 
presenting the computations for the curvature evolution equations. 



5.1 A simple example 

In this example we will suppose that we have a vector that evolves along 
the world line of the central vertex according to 

— - NF^ (5.1) 
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Our aim is to obtain a related equation that describes the evolution of the 
vector along the whole length of the world line, not just the short section 
contained within this one cell. 

Suppose that we have an orthonormal tetrad Ca — e^ad^i a — 1,2,3,4 on 
u with ei aligned to the future pointing normal to u, and that we 

have aligned the RNC coordinate axes with the tetrad (note how this gives 
precedence to the tetrad over the coordinates) . Thus at the central vertex of 
Q we have 

da . C^a ^^a i ^-^u ^^i 

Qi^^ = diag(-l, 1, 1, 1) , Qab = diag(-l, 1, 1, 1) 

We now propose the following evolution equations along the world line of the 
central vertex in Q. 

^ = e^VW, ^ = -e.V.^ (5.2) 

where V^A^ = {±N^^)e''i and VW = {±N^'')e„\ i = 2,3,4. What can we 
say about the evolved data? First, note that the orthonormal conditions are 
preserved, that is 

dt' ~ ' dt' 

Thus the tetrad obtained by integrating the above equations will remain 
orthonormal along the world line of the central vertex. Second, using 

(Nn^").^ = N^n^" - ±(N''')n^ - NK^"^ (5.4) 

to compute dn'^/dt' — n'^.y {Nn") we see that 



de'^i dni^ de// dn^ 



(5.5) 



which shows that e'^i = and e^^ = — everywhere along the world line. 
That is, e^i remains tied to the world line. All that remains is to account for 
how the tetrad rotates around the world line. This we shall do by evolving the 
projections of the e'^j, i = 2, 3, 4 onto the legs of the lattice. Let Va = v^ad^, 
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a = 1,2,3 be any three distinct legs of the lattice attached to the central 
vertex. Now consider a short time step in which the vector Va sweeps out a 
short quadrilateral in spacetime (see figure (2)). The upper and lower edges 
will be the past and future versions of Va while the remaining two sides will 
be generated by the word lines of the vertices that define Va- Since we have 
assumed at the outset that all vertices evolve normal to the Cauchy surface 
we see that these vertical vectors correspond to Nn^. The important point 
is that this set of four vectors forms a closed loop, in short the vectors Va 
and Nn^d^ commute, thus 

v''a;u{Nn^)=v\{Nn^').^^ (5.6) 

The left hand side is simply dv^/dt', while the right hand side can be ex- 
panded using (5.4). This leads to 

^ = {N^n^-NK^,)v\ (5.7) 

where we have dropped the term involving n^v^a as this would be O (L"^) 
with m > 2 while the remaining terms are all O (L). 

We are now ready to construct our scalar evolution equations. Let Wa := 
IV^e^a and Va'' := v^'ae^,^ then 



dt ~ dt' ^ dt' 



^ = ^^e/ + w'^a— ^ , 6 = 1,2,3,4 
dt' dt' ^ dt' ' ' ' ' 

Each of these equations can be re-cast entirely in terms of the scalars by first 
using (5.2,5.3,5.7) to eliminate the time derivatives on the right hand side 
followed by the substitutions = Wae^"" and v^a = vj'e^b- This leads to 



dt 

dW,, 



dt' 



NTi + WnViN, 2 = 2,3,4 (5.9) 
dvj IdN , 

- jVa (5.10) 



dt' N dt 

dv, 



'a 



dt' 



NV^vJ , t, J = 2,3,4, a = 1,2, 3 (5.11) 
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where we have introduced the scalars W„ = VF^n'^, Tn = F^n^, Ti = F^e^i, 
and lOj = K^yC^^e'^ j. These are our final equations. They are vahd along 
the whole length of the world line, not just the part contained in one cell. 

Equation (5.11) describes the motion of the tetrad relative to the legs of the 
lattice. As we integrate forward in time we can use the values of Va to locate 
the tetrad within the computational cell. If we chose to construct an RNC 
within the cell then we can go one step further and recover the values of e^i 
and the W^^. 



5.2 Curvature evolution equations 

Now we can return to the task of constructing the generalised evolution 
equations for the curvatures. We start by introducing a pair of relations 
between the tetrad and coordinate components of the curvature tensor 

"T^abcd RfiavP^^ d 

and then forming a typical evolution equation 

dT^abcd dR^gyfs a d {e^ae^be\e^d) 

= ,e ,e ,e . + (5.12) 

with each d/dt' term on the right hand side replaced by a suitable combina- 
tion of the existing evolution equations, (4.4-4.17) for the curvature terms 
and (5.2,5.3) for the tetrad terms. 

Rather than working through all 14 equations we will demonstrate the pro- 
cedure on just one equation (4.4) leaving the remaining equations (but not 
their working) to the Appendix. So our starting point is 

dV dJ^ d (p^^ p" 

^i^xyxy 'J'-n.xyxy , jd \^ x^ y^ x^ y ) 

~dF~ ^ + dt' 

and using (4.4) we obtain 

dV d (pi^ p^ 

U'l^xyxy J-, J-, J-, 1^ xC^ yC^ xC^ yj 



Rtyxy,x Rtxxy,y ~l~ RfiavP 



Finally we use (5.2,5.3) to eliminate the time derivative of e^a, leading to 
dn. 



"xyxy 



~ Rtyxy,x Rtxxy,y 

+ T^tyxy^ xN — TZtxxy^ yN + TZtyxy^ xN — TZtxxy^ yN (5.13) 
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This is as far as we need go, though it is tempting to make the substitutions 
Rtyxy = 'Tlabcdet"'ey^e^''ey'^ and Rt^xy = 'Tlabcdet'eJ'e^''ey'^. But that is not 
really necessary as we can defer those substitutions until we actually need 
values for the stated partial derivatives. This is described in more detail in 
section (7). 

Note that when introducing the lapse function by the substitution t = Nt' 
we have not made explicit the coordinate transformation on the curvatures 
other than to use distinct labels t and t'. In this way we use t' as an integra- 
tion parameter on the world line of each vertex while retaining the original 
coordinates {t, x, y, z) as the local Riemann normal coordinates (and thus at 
any point on the world line we continue to have {g^u)o = diag(— 1, 1, 1, 1)). 
We choose to maintain this distinction between t and t' not only to keep the 
equations tidy but also because it leaves the equations in a simple form well 
suited to numerical integrations. 

Clearly the above procedure can be applied directly to each of the remaining 
13 curvature evolution equations. The final results for all 14 equations can 
be found in the Appendix. 

5.3 Hyperbolicity and constraint preservation 

It is natural to ask if the new system of evolution equations are hyperbolic 
and also, are the new constraints preserved by the new evolution equations? 
The answer to both questions is yes and we will demonstrate this as follows. 

Given that IZabcd = R^lav|3^^ a^°'bG^ c^^ d we see that 

7^afecd,e/ = R^.a.p,pre\e''be\e^de'ee^f + V^bcdef {R, N, dR, dN, d'^N) 

where Vabcdef is a function of R^aui3, N and the indicated partial deriva- 
tives. Importantly, Vabcdef does not contain any second partial derivatives 
of the curvatures. We have previously shown that, at the central vertex, 
each R^avp satisfies a wave equation of the form = g'^'^ R^aui3,pT with g^'^ = 
diag(— 1, 1, 1, 1). Thus we find that 

g'^nabcd,ef = g'^Vabcdef {R, dR, N, dN, d^N) 

where g'^-^ = diag(— 1, 1, 1, 1). It follows that each TZabcd satisfies a wave equa- 
tion with source terms and therefore we have shown that the new evolution 
equations constitute a hyperbolic system. 
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A similar analysis can be applied to the constraints. We begin by writing a 
typical differential constraint (4.28-4.33) in the form 

= W^^,p{dR) 

where the right hand side depends only on the the first derivatives of R^avp- 
Introducing the lapse function is trivial (there are no time derivatives, so the 
equation is unchanged). If we define the frame components Wabcd by 

then we find 

and as we have previously shown that W^aup = and W^aui3,p = it follows 
that Wabcd = and Wahcd,t = 0. It is easy to see that the same procedure can 
be applied to the remaining constraints (4.18-4.24) with the same outcome. 
Thus we have shown that the new constraints are conserved by the new 
evolution equations. 

6 Coordinates 

There are at least two instances where the vertex coordinates are required. 
First, when constructing the transformation matrix used when importing 
data from neighbouring cells. Second, as part of the time integration of leg- 
lengths, equations (3.1-3.2). They are also required when computing the 
extrinsic curvatures (7.1) and the hessian (7.2). 

Recall that within each cell we employ two distinct coordinate frames, one 
is tied to the tetrad associated with the central vertex while the other is 
aligned with the lattice. Both frames share the central vertex as the origin. 
We will describe first how to construct the lattice coordinates, which we will 
denote by y^, followed by the tetrad coordinates, denoted by x^. The lattice 
coordinates are only ever used in the construction of the tetrad coordinates, 
once these are known then the lattice coordinates can be discarded. Note 
that terms such as Rxyxy, K^y^z etc. are referred to the tetrad coordinates. 

For a large part of this discussion we will be concerned mainly with the 
scaling of the coordinates with respect to the typical lattice scale (e.g., to 
establish that t = 0{L^)). This apphes equally well to both coordinate 
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frames and so, to be specific, we will present the arguments in terms of the 
tetrad coordinates. Once we have sorted out these scaling issues we will 
compute the lattice coordinates followed by the tetrad coordinates. 

Our first task will be to construct the piece of the Cauchy surface that is 
covered by a typical computational cell. Recall that we view the Cauchy sur- 
face to be a smooth 3-dimensional surface that passes through each vertex 
of the lattice and that it shares with the lattice, at each vertex, the same 
future pointing unit normal and second fundamental form (the extrinsic cur- 
vatures). In our local Riemann normal coordinates we wish to construct an 
equation of the form = —t + f{x^) that passes through the vertices of this 
computational cell and with given extrinsic curvature at the central vertex. 
For this we use the familiar definition that Sn'^ = —K^^Sx^ for the small 
change in the unit normal under a displacement across the Cauchy surface. 
If we take the displacement to be from the central vertex (o) to a nearby 
vertex (a) then we have 

<-< = -ir^< (6.1) 

But we chose the coordinates so that = (1, 0, 0, 0)^ while for the surface 
= — t + f{x^) the unit normal at (a) is simply n^' = g'^'^{—l,f^u)i,/M = 
(1, f^uY/M where M = 1 + O (L^) is a normalization factor. Thus we have 
(1, f jy = (1, 0, 0, O)'' - K^^x"^ + O (L2) and this is easily integrated to give 

ta = -\K,,x^a< + 0{L') (6.2) 

Note that since K'^yU^ = we can use this last equation to compute the 
time coordinates for each vertex in the computational cell (given the spatial 
coordinates x^ and the extrinsic curvatures K^v)- 

Consider the geodesic segment that joins the central vertex (o) to a typical 
nearby vertex (a). Then from the definition of Riemann normal coordinates 
we have 

< = m^aLoa (6.3) 

where is the unit tangent vector to the geodesic at (o).^ Thus it follows 
that 

\<\ = 0{L) (6.4) 

for each vertex in the computational cell. Combining this with the above 
equation (6.2) for ta shows that 

\ta\ = 0{L') (6.5) 

^Actually, by virtue of the fact that the path is a geodesic segment expressed in Riemann 
normal coordinates, the values for m'^ are constant along the geodesic. 
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This result could also be inferred from the simple observation that m* — )■ 
as L — )■ (this is a consequence of the smoothness of the Cauchy surface at 
(o)). 

We turn now to the simple question - How accurate do we need the coor- 
dinates to be? That is, if are the exact Riemann normal coordinates for 
vertex i, then how large can we allow — | to be? The answer can be 
found by a simple inspection of the evolution equations (3.1-3.2). The trun- 
cation terms in those equations are O (L^) thus we can safely get by with 
O (L^) errors in the coordinates, that is 

|<-xf| = 0(L^) (6.6) 

The good news is that such coordinates are readily available - flat space will 
do. To see that this is so, assume, for the moment, that we have estimates 
for the ii'^t/ and then look back at equations (6.2,1.2). This is a coupled 
system of equations for the coordinates (t,x,y,z)'^ for each vertex in the 
computational cell. We are fortunate to have an explicit equation for the 
time coordinates, namely (6.2). This allows us, in principle, to eliminate 
each time coordinate that appears in equation (1.2). The result would be a 
set of equations for the spatial coordinates x". In the following we will not 
make this elimination explicit but take it as understood that such a process 
has been applied. We will have a little more to say on this matter in a short 
while. 

For a typical vertex (/) we will need to compute three spatial coordinates and 
thus we look to the legs of a tetrahedron. Suppose that that tetrahedron has 
vertices (ijkl) and suppose that we have computed, by some means, the exact 
Riemann normal coordinates x'^ for vertices (ijk). The exact coordinates xf 
for vertex (/) could be obtained by solving the system of equations 

but we could also construct flat space coordinates for vertex I by solving 
the system 

Ll = 9,AS^'a - <Wa - A) a = k (6.8) 
From the last equation we conclude that — = O iV) for a ^ I. Next, 
make the trivial substitution xf = xf + (x^ — xf) in the flrst term in (6.7), 
expand and use (6.8) to obtain 

= —2g^y{x'^ — xf )(x^ — x\) + gfi,y{xf — xf){x'j^ — x\) 

— —Rfj^ai^iijx'^x^Xi xf a = i,j, k 



19 



and as each = C (L) for a = k we easily see that 

\xr-xr\ = 0{L') (6.9) 

The fly in the ointment in the above analysis is the assumption that we 
knew the i^^j^ (and thus we could eliminate the ta). This is not exactly 
correct for the K/j^u are found by solving equations (3.1) which in turn requires 
the coordinates x]^ which we have yet to compute (at that stage). Luckily, 
this is not a major problem. Look carefully at equation (6.7) and recall 
that g^i, = diag(— 1, 1, 1, 1). Thus the t-terms will appear only in the form 
— {ta — tiY and in the curvature terms of the form RtuvwtaX^x^xf . The point 
to note is that since t = O (L^) we see that each of these terms is O (L") 
with n > A and thus they have no effect on the above analysis. Thus even 
though we argued previously that we should eliminate the ta using equation 
(6.2) the above argument shows that we can put ta = without harm. 

Our final calculation concerns the errors induced in ta by using the approxi- 
mate Xa and Kfj^iy rather than their exact counterparts. Our analysis is very 
similar to that just presented. We start with the two sets of equations, the 
approximate and exact equations, 

2ta = -Ku.xyi and 2^ = -iT^xX (6-10) 

We will assume that \Kuv — Kuv\ is at least O (L) (this is one assumption 
that we will not relax at a later stage). Then we make the trivial substitution 

xf = xf + (xj* — x") as above to obtain 

^Ky^y Ky^J^ XX 2K yyX a {S^a •^a) -^uv iS^a ■^a) i.^ a "^a) 

(6.11) 

Using xl = 0{L),S:l = (L) and - Kyy\ = O (L) we find that 

\ta-ta\ = 0{L'') (6.12) 

6.1 The lattice coordinates 

We return now to the concrete question of how to compute the vertex co- 
ordinates within one computational cell. We will first compute the lattice 
coordinates followed by the tetrad coordinates x^. Our present challenge 
is to find the solutions of the coupled system of equations 

Ll, = g,.{y^a-yWa-yl) (6.13) 



20 



for a suitable subset of the legs (ab) in the computational cell (equal in num- 
ber to the number of unknown coordinates). The problem here is that if we 
treat this as a system of equations for the spacetime coordinates {t, x, y, z)'^ 
it is extremely unlikely that we will find any solutions (or if we do then the 
numerics will almost certainly be extremely unstable). The reason is quite 
simple - the vertices are assumed to lie within one 3-dimensional Cauchy 
surface. This suggest that we should use the above equations to determine 
the spatial coordinates {x,y,z)^ with the time coordinates found by other 
considerations. Fortunately we already know, from the above analysis, that 
each \ta\ = O (L^) while \y'^\ = O (L). Thus we see that all terms involving 
the ta are O (L^) and thus will be consumed by the O {L^) truncation errors 
inherent in the above equation (as an approximation to equation (6.7)). So 
we may safely discard all the of the ta terms in the above equations. The 
next trick that we will use is the observation that the coordinates can be 
computed one vertex at a time. This is easily shown by direct construction. 
Consider a typical tetrahedron with vertices (oijk) where (o) is the central 
vertex and suppose we have computed the coordinates for (ojk). Our task 
now is to solve the following equations 



where the last pair of equations were obtained by expanding L^^ = QuviVa " 
yb)hJa ~ Vb)- ^ simple calculation shows that the solution is given by [4] 




guvVkUk 

'^guvyM 



(6.14) 
(6.15) 
(6.16) 



yl = Py't + Qy] + Rn 



u 



where 



,uv xyz r s 



P 




mjkrriij 



Q 




(t2 _ p2r2 



Q^L% - 2PQm,,) 



1/2 



R 



± 




oi oj 





2mik 




2mjk — LIj + LIi^ — L?^^ 
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The two solutions, one for each choice of the ± sign, correspond to the 

two possible locations of the third vertex (A;), one on each side of the plane 
containing the triangle (oij). Which choice is taken will depend on the design 
of the lattice. A systematic choice can be made by noting that the vectors 
l/f , and form a right handed system. With R > the vector lives 
on the same side of the plane as n". 

To complete the picture we need coordinates for the first two vertices (1) 
and (2). Since we chose to align our coordinates so that the x-axis passed 
through vertex (1) while the vertex (2) is contained in the xy-plane we must 
have Ui = {A, 0, 0)" and ^2 — {B, C, 0) for some numbers A > B and 
C > such that 

-^01 — QuvViyi 

-^02 ~ 9uvy2y2 
-^01 + -^02 ~ -^12 ~ 2(7u„|/"|/2 

The solution is readily found to be ^4 = Lqi, B — [L"^^ + — L\^)/{2Lqi) 
and C = - B^y/\ 



6.2 The tetrad coordinates 

The transformation from the lattice to tetrad coordinates is quite simple. 
Let Ca be the basis for the tetrad frame and let be the corresponding 
basis for the lattice frame. Recall that we have previously chosen the frames 
so that both ei and dt are aligned with the normal to the Cauchy surface. 
Now consider a typical vector Va that joins (0) to (a). In the lattice frame 
this vector has components while in the tetrad frame, with basis e^, its 
components are just Vg,^. That is we have, for a = 1, 2, 3 

n^dt^ei (6.17) 
Vl = vj (6.18) 

Va = y^d^ = Va'cb (6.19) 

In the last equation both the and vj' are known. Thus we have sufficient 
information to compute (9^ in terms of Cq and vice versa. Note that the tetrad 
coordinates are given by 

x'^^Va'e'^b with e\^S^ (6.20) 
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Finally, using equation (6.2), we can compute the time coordinate for every 
vertex, not just the three vertices associated with Va ^ a = 1, 2, 3 



a = 1, 2, 3, ■ ■ ■ 



(6.21) 



7 Source terms 

We have previously mentioned, without giving details, that source terms such 
as Rxyxy,z can be computed by applying a finite difference approximation to 
data imported from neighbouring cells. Here we will outline how such a 
procedure can be applied (the exact details will of course depend on the 
structure of the lattice). The same procedure can also be used to estimate 
the spatial derivatives of the e^a- 

Suppose we have two neighbouring computational cells that have a non-trivial 
overlap (as indicated in Figure (1)). Each cell will carry values for Rxyxy in 
their own local RNC frames. Our first task would be to import the values 
form the one cell to the other. This will entail a coordinate transformation, 
composed of a boost (to account for the change in the unit normal between 
the two cells) and a spatial rotation (to account for the different orientations 
of the legs of the cells). 

Let be the (tetrad) coordinates in one cell and let x'^ be coordinates in the 
other cell. Our plan is to import data form the x'^ frame to the x^ frame. We 
will demand that the overlap region be such that it contains at least one set 
of three linearly independent vectors (i.e., legs), at O', which we will denote 
by Wj, i = 1,2, 3. Since we know the coordinates of each vertex in each cell 
we can easily compute the components of Wi, z = 1, 2, 3 in each frame. The 
normal vector Uo' at O' will have components n'^ = (1,0,0,0)^ in the x''^ 
frame. But in the x^ frame we expect n'^, = ra^ — K'^^.x'^,. Thus we have 4 
linearly independent vectors at O', expressed in two different frames, and so 
there must exist a mapping from the components in one frame to those in 
the other. That is there exists a If^i, such that 



Since we have values for the components of Uo' and Wi, i = 1,2,3 in both 
frames we can treat this as a system of equations for the If^u- 



< = W'.w'^ , z = 1,2,3 



(7.1) 
(7.2) 
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With the U^,y in hand, we can compute the values of R^aui3 at O' in the x'^ 
frame of O by way of 



(7.3) 



with Uf^" = gfj.aQ^^U'^ p and g^y = diag(— 1, 1, 1, 1). This can be repeated for 
all of the vertices that surround O. The result is a set of point estimates for 
R^avjs in the neighbourhood of O which in turn can be used to estimate the 
derivatives of R^aup at O. This part of the process is similar to that required 
when computing the Hessian (see below) and presumably similar methods 
could be applied. 

Note that for a sufficiently refined lattice, the f/^j^ should be close to the 
identity map, that is U^i, = 6^^, + V^^O (L) where the V^i, are each of order 
O (1). This can be used to simplify some of the above computations. 

See [1] for a complete example in the context of the Schwarzschild spacetime. 

In section (5.2) we noted that substitutions such as Rtyxy = Tlabcd^t^^y^^x'^^y'^ 
could be introduced into the curvature evolution equation (5.13). At that 
time we argued that that was not necessary for the coordinate data, in this 
instance Rtyxy ^ could easily be recovered when needed by using Rtyxy = 
T^abcdet°'ey''ex'^ey'^ . Then the scheme described above could be used to com- 
pute Rtyxy, X- However there may be numerical advantages in making a formal 
substitution before estimating any of the partial derivatives. For Rtyxy,x this 
would lead to the following 



Since the TZabcd are scalars, their partial derivatives can be estimated without 
requiring any of the frame transformations described above (importing such 
data from neighbouring cells is trivial). This leaves us with the derivatives 
of the form (e^")^^;. Since = —e^^ we can use (5.4) to eliminate any of 
the spatial derivatives of e^^, in this case (e^^)^^,. This would introduce the 
extrinsic curvatures into the evolution equations. However the remaining 
partial derivatives, (e^*)^^;, i = 2,3,4, would have to be estimated using the 
methods described above (by importing data from neighbouring cells etc.). 
This approach does incur a small computational overhead which may be 
justified if it brings some improvement to the quality of the numerical data 
(e.g., better accuracy and or stability). Judging the merits of this variation 
against the simple method given in section (5.2) might best be decided by 
direct numerical experimentation. 



R, 



'tyxy,x 



ij^abcd^t ^y ^x 



■abcd,x^t ^y ^: 
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7.1 Extrinsic curvatures 



A cursory glance at equation (3.1) might give the impression that it consti- 
tutes a simple linear system for the Kuv But things are never as simple as 
they seem. The problem, as already noted, is that there are far too many 
equations for the six K^v If we make the reasonable assumption that the 
lattice data is a good approximation to the (unknown) continuum spacetime 
then we can expect considerable redundancy in this overdetermined system. 
How then do we pull out just six equations for the six Kuv"^ One option is to 
reject all but six of the equations and hope that this yields an invertible sys- 
tem for the Kuv A better, and more flexible approach, is to take a weighted 
sum of the equations, that is we create a new set of equations of the form 

= ^P-^ - Ku.Ax:,Ax:,) (7.4) 

ab 

where VF"^ are a set of weights of our own choosing (typical values being and 
±1). With n = 1, 2, 3 ... 6 we have six equations for the six unknowns. This 
idea has been used previously [4] and worked very well. There are certainly 
other options that could be explored (e.g., different choices of weights, least 
squares) but we have tested none simply because the above scheme seems to 
work well. 



7.2 The Hessian 

At some point we will need to estimate the N\uv at a central vertex. Since 
iV is a scalar function and since we are using Riemann normal coordinates 
this computation is essentially that of computing all of the second partial 
derivatives on an unstructured grid. There is an extensive literature on this 
point in the context of finite element schemes. We mention here one approach 
which we discussed in one of our earlier papers [ I ] (but which we have yet to 
test). 

Consider a typical leg {ij) in some computational cell. We can estimate N\u 
at the centre of the leg by the centred finite difference approximation 

in which {mu)^^ is the unit vector tangent to the geodesic and oriented so 
that it points from (i) to (j). We can repeat this computation for each leg 
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in the computational cell and then estimate N\u^ by a least squares fit of the 
function 

iV|„(x) =iV|, + iV|„,a;'' (7.6) 

to the data generated above by equation (7.5). A suitable least squares sum 
would be 

u ij 

where x^j is the centre of the leg (ij). Note that this least squares fit must 
be made subject to the constraint N\uv = N\^u- The coefficients N\u and N\uv 
would then be taken as our estimates for the corresponding quantities at the 
central vertex. 



8 Discussion 



There are a number of aspects of this paper that could easily be debated. 
For example, should we proceed with the substitutions such as Rtyxy = 
'T^abcd^t'^y'^x^^y'' in equatiou (5.13)? As already noted in section (7) this 
would introduce a raft of new terms including the extrinsic curvatures. We 
chose not to use the substitution solely for reasons of simplicity. There is also 
a question over our choice of tetrad. Do we really need to demand that the 
tetrad be orthonormal? Not at all. We could choose to tie the tetrad to the 
legs of the lattice (and then the tetrad would no longer be needed) but that 
would produce a coupling amongst all of the evolution equations (e.g., the 
evolution equation for TZxyxy would be a linear combination of all of the evo- 
lution equations for R^avp)- The resulting equations would not be anywhere 
near as simple as those listed in the Appendix. Then we have the issue of 
estimating partial derivatives on an irregular lattice (for the Hessian and the 
source terms in the curvature evolution equations). This is non-trivial but at 
least there is an extensive literature on the subject and so a workable solu- 
tion should not be too hard to find (which may be the least squares method 
suggested in section (7.2)). All of these issues (and most likely others) can 
be explored by direct numerical exploration on a non-trivial 3 -f 1 spacetime. 
We plan to report on such investigations soon. For a simple application to 
the 1 + 1 Schwarzschild spacetime see [1]. 
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A The curvature evolution equations 



Here we list all 14 curvature evolution equations (this follows on from sec- 
tion (5.2) where we provided details of the derivation for the first equation 
below). 



dT^xyxz 

dt' 



dT^xyyz 



dT^xzxz 

dt' 
dTZ 



dT^yzyz 



dU 



■tyxy 



dT^txxy 



dT^tzxy 



dn 



tzxz 



dt' 



Rtyxy,x Rtxxy,y 

+ TZty^yV^N - UtxxyVyN + UtyxyV^N " UtxxyVyN (A.l) 

Rtzxy,x Rtxxy,z 

+ Htyxz^^N - TZtxxzVyN + TZtzxyV - TZtxxyV.N (A.2) 



^tzxy,y Rtyxy,z 

+ 'R-tyyZ^xN - Tlt^yzVyN + TZu^yVyN - TZty^yV.N (A.3) 



Rtzxz,x Rtxxz,z 

+ 7^^.x.V,A^ - TltxxzV.N + TltzxzV^N - UtxxzV ,N (A.4) 



xzyz j-y p 



J^tzxz,y J^tyxz,z 
+ TltzyzV^N - TltxyzV.N + UtzxzV yN - UtyxzV ,N (A.5) 



^tzyz,y Rtyyz,z 

+ TltzyzVyN - -Rtyy^V-^N + TZt^y-NyN - TZtyy.V.N (A.6) 



+ TZiy^yV'N + TZtytyV^N - TZtxtyVyN (A. 7) 



Rxvxv.v Rx 



^^xyxy,y J-^xyxz,z 
+ Ui^^yV'N + UtxtyV^N - UtxtxVyN (A.8) 



RxAIXZ.X -^.7 



^xyxz,x ~r J^xyyz,y 

+ ni,,yV'N + -RtytzV^N - TZt^t^VyN (A.9) 



Rxzxz,x ~l~ Rxzyz,y 
~l~ T^izxz 

V'N + nt,tzV,N-nt,tzV,N (A.IO) 
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dlZtxxz 



dT^tzyz 



dT^tyyz 

dt' 



Rxyxz,y Rxz 



^^1 ^^xyxz,y 
~l~ T^ixxz 

V'N + ntxtzVxN -UtxtxV.N (A.ll) 



+ TZiy^.V'N + -RtytzVxN - TZtxtyV.N (A. 12) 



Rxzvz.x ~l~ Ry 



~ ^^xzyz,x ^ J'^-yzyz,y 

+ TZi.y.V'N + TZt^tzVyN - TZtytzV.N (A. 13) 



RxyyZyX Ryzyz,z 



+ Tliyy^VN + TltytzVyN - -RtytyVzN (A. 14) 



Note that in the above there are two instances of TZtxyzi in (A. 3) and (A. 5), 
and these should be replaced with IZtyxz — T^tzxy 



B Riemann normal coordinates 



We recall here a few basic properties of Riemann normal coordinates. A set 
of coordinates are said to be in Riemann normal form if every geodesic 
passing through a given point O (the origin) is described by x'^(s) = sv'^ 
where s is an affine parameter and f'^ is constant along the geodesic. It 
follows from the geodesic equation and its successive derivatives, that the 
connection and its higher symmetric derivatives" all vanish at the chosen 
point, that is at O 

= r:^^., (B.i) 

= rf„,„,;.3...«„) n = 3,4,5,--- (B.2) 

These conditions do not uniquely determine the coordinates for we are free 
to apply a transformation of the form ^ A'^^x'^ which clearly preserves 
the property that the geodesies through O are of the form x'^{s) = sv^. 
This freedom can be used to ensure that the metric at O is simply g^i, = 
diag(-l, 1,1,1). 

Choosing the coordinates so that the connection vanishes at the origin does 
introduce some nice properties, in particular covariant differentiation reduces. 



"Here we take a small liberty with notation, the upper index on the Christoffel symbol 
should be ignored when computing covariant derivatives. 
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at the origin, to simple partial differentiation. This fact was essential to the 
analysis given in sections (4). 

There are two main impediments to the existence of Riemann normal coordi- 
nates. The metric must be smooth throughout the neighbourhood (i.e., away 
from curvature singularities) and each point in the neighbourhood should be 
connected to the origin by exactly one geodesic (i.e., no pair of geodesies 
through O should cross, except at O). These conditions are easily satisfied 
by simply choosing the neighbourhood around O to be sufficiently small (but 
not vanishingly small). 

In these coordinates the metric and connection can be expanded as a Taylor 
series around O leading to 

g^Lu{x) = Qf^u - ^R^iaupX'^x'^ - ^R^auPn^'^x'^x'^ + O (L^) (B.3) 

g^^-[x) = g'^" + ^-R^J'px'^x^ + ^-R^^^p^.x'^x^x'^ + O (L^) (B.4) 

+ {a^ P) + [L^) (B.5) 

If we know the Riemann normal coordinates, x^ and x^, for a pair of points, 
i and j, then we can compute the length of the geodesic segment that joins 
the points by 

^Ij = (^Qfiu — -^R^J.au|3x1jX^j — -R^aui3,'rX"jX^jX]j^ Ax'^jAx^j + O (L^) (B.6) 

where Ax^- := x^ — xf and x^- := (x^ + xf)/2 is the mid-point of the leg. 
The unit tangent vector m^- to the geodesic at i, is given by 

L.^mf^. = Axf^. + ^x°Axr, Axji?^„;3 + ^x-x'^AxjAxZ^.i?^.,,^ 

+ ^x"x'^Axf,.Ax7,.i?'^;3.,,„ + ^x-x'^AxgAxZ^.i?^^^ (B.7) 

+ ^x-Axr^.AxgAxXi?^;,,,,, 

Finally, if we have a geodesic triangle built on the three points z, j, k then 
the generalised cosine law takes the form 

2L,fcL,fcCos% = Ll+L%-Ll-^-R^^,pAx>:,Ax^,Ax^,Ax% + {L') (B.8) 

in which 9ij is the angle subtended at vertex k by the geodesic that connects 
i to j. 
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Figure 1: An example of the overlap, the shaded region, between a pair of com- 
putational cells. The central vertex of each computational cell is denoted by the 
large dots whereas the smaller dotes denote the vertices that define the bound- 
ary of the computation cells. These vertices are themselves the central vertices of 
other computational cells. In this 2-dimensional example the overlap consists of 
just the pair of triangles. In 3 dimensions the over lap would consist of a closed 
loop of tetrahedra. In each case there is ample information available to obtain a 
coordinate transformation between the pair of local Riemann normal frames. 
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{Nn)oSt' 



O 



V 



a 



{Nn)a5t' 



a 



Figure 2: Here we show the evolution of one leg (oa) within one computational 
cell. Clearly the four vectors form a closed loop and thus {Nn)o5t' + v'^ = Va + 
{Nn)a)6t' which leads directly to equation (5.6). 
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